/* Poisson Equation solver using JACOBI iteration. */
# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <stdlib.h>

# ifdef _OMP_SOLVE_POISSON_OMP_
# include <omp.h>
# endif
/*
# include <time.h>
# include <sys/timeb.h>*/

/* This structure is used to detect most sensitive cell in an itteration. */
struct dpmax
    {
     short im, jm, km;
     double dpm;
    };

extern void press_bc(const struct grid *g, short tflag, unsigned short iRK);

extern long int solve_poisson_omp(struct grid *g, const double ***d, unsigned short iRK, double ***fspace)
  {
   int cnt;
   short i, j, k, nr, nth, nz, stop_cnt, tid, im, jm, km;
   double dp, ***p, *r, *th, *z, prf, ***pmemloc, dpm;
   char flag;
   struct dpmax *tdp;

//  struct timeb t1,t2;

#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'solve_poisson_omp()'...\n");
#  endif

   r = g->r;   z = g->z;   th = g->th;   nr = g->nr;   nth = g->nth;   nz = g->nz;

   prf = PRF_INIT; stop_cnt = 0;
   tdp = (struct dpmax *)malloc(sizeof(struct dpmax)*g->nCPU);

   pmemloc = g->p;
   flag = 'o';
// ftime(&t1);
   for(cnt = 1; cnt <= POISS_CNT_MAX; cnt++)
      {
       if(flag == 'o'){ p = pmemloc; g->p = fspace; flag = 'e';}
       else{ p = fspace; g->p = pmemloc; flag = 'o';}

       tid = 0;
    /* Iterating through all the cells using JACOBI method. */
#      ifdef _OMP_SOLVE_POISSON_OMP_
       #pragma omp parallel private(i,j,k,dp,tid)
#      endif
           {
         /* Initializing elements of the structure dpmax. */
#           ifdef _OMP_SOLVE_POISSON_OMP_
            tid = omp_get_thread_num();
#           endif
            tdp[tid].dpm = (double)0;
            tdp[tid].im = -1;
            tdp[tid].jm = -1;
            tdp[tid].km = -1;

#           ifdef _OMP_SOLVE_POISSON_OMP_
            #pragma omp for schedule(guided)
#           endif
            for(k = NzG; k < NzG+nz; k++) for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
                 {
                  dp = P_ijk_NEW - p[i][j][k];

                  if((dp<0?-dp:dp) > (tdp[tid].dpm<0?-tdp[tid].dpm:tdp[tid].dpm))
                      {
                       tdp[tid].dpm = dp;
                       tdp[tid].im = i;
                       tdp[tid].jm = j;
                       tdp[tid].km = k;
                      }
                  g->p[i][j][k] = p[i][j][k] + prf*dp;
                 }
           }

    /* Calculating the absolute maximum. */
       dpm = tdp[0].dpm;
       tid = 0;

       for(i = 1; i < g->nCPU; i++)
         if( (tdp[i].dpm<0 ? -tdp[i].dpm : tdp[i].dpm) > (dpm < 0 ? -dpm:dpm) )
             {
              dpm = tdp[i].dpm;
              tid = i;
             }
       im = tdp[tid].im;
       jm = tdp[tid].jm;
       km = tdp[tid].km;

       press_bc(g, 0, iRK);

       if(cnt == 1 || cnt%DISP_FREQ == 0)
            {
             fprintf(g->dptr, "=> %d: dpm = %E(%d, %d, %d)\n", cnt, dpm, im, jm, km);
             fflush(g->dptr);
            }

       if( (dpm < 0 ? -dpm:dpm) <= PRESS_DELTA) ++stop_cnt;
       else stop_cnt = 0;

       if(stop_cnt == POISS_CNT_MIN) break;
      }

   free(tdp);

/* If the new pressure field is located at "fspace", copy it to the memory located at _p. */
   if(flag == 'e')
      {
       for(i = 0; i < 2*NrG+nr; i++) for(j = 0; j < 2*NthG+nth; j++) for(k = 0; k < 2*NzG+nz; k++)
       pmemloc[i][j][k] = fspace[i][j][k];

       g->p = pmemloc;
      }

   flag = (stop_cnt == POISS_CNT_MIN ? 'Y':'N');

   fprintf(g->dptr, "N_ittr = %d; dpm = %E(%d, %d, %d); Conv. = %c\n", cnt, dpm, im, jm, km, flag);
   fflush(g->dptr);

   fprintf(g->tptr, "N_ittr = %d; dpm = %E(%d, %d, %d); Conv. = %c; ", cnt, dpm, im, jm, km, flag);
   fflush(g->tptr);

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'solve_poisson_omp()' ends.\n");
#  endif

   return cnt;
  }
